Active neurons produce small currents whose basis is the change in the concentration of certain ions (ionic currents). The flow of current in the neural system produces a weak magnetic field. The measurement of this field outside the head and the estimation of the neuronal current which produced this field is called Magnetoencephalography (MEG). In spite of the fact that neuromagnetic signals are only of the order of 50-500 fT, it is still possible to measure them using the superconductivity quantum interference device. Electroencephalography (EEG) measures the electric potential on the scalp. Although EEG plays a significant role in clinical practice especially in epilepsy, it is not routinely used to produce images depicting the areas of the brain which are more active; such images are produced via MEG. In order to produce images of brain activation using either MEG or EEG it is necessary to solve certain mathematical inverse problems. The relevant inverse problems for MEG and EEG involve the calculation of the neuronal current from the knowledge of the magnetic field and of the electric potential, respectively. Currently there exist two main approaches for the solution of the above inverse problems: 1. Dipole Models. Here it is assumed that the neuronal current is localized only at a small number of points. At each of these points the current is expressed in terms of a dipole characterized by the location Tj and the moment Qj, j=1; 2; . . . N, where N is small, and the constants {Qj,τj}1N are to be determined. The magnetic field and the electric potential are computed starting with the above current, and then the constants {Qj,τj}1N are adjusted step by step minimizing a certain natural cost function. 2. Distributed Models. Here the continuous current is approximated again via a sum of dipoles, but now N is very large of the order of 1000 to 10,000. Then, the procedure is similar to the one described for dipole models.
The first, dipole models approach described above is based on the assumption that a small number of dipoles is sufficient for modelling the neuronal current. The unknown parameters are the positions and orientations of these dipoles and these are adjusted to minimize the mismatch between the values of the electric potential in the case of EEG, and of the magnetic field in the case of MEG generated by the dipoles and the actual measurements. However these methods rely on local optimisation and suffer from getting trapped by local minima. Furthermore, they require the assumption that the number of dipoles is less than the number of sensors. The second, distributed source modelling approach divides the cerebrum into voxels each representing a possible of location of a current source (there is no a priori assumption regarding the number of dipoles), and the objective is to find a configuration of activity at these locations which is consistent with the measurements. However there exists an infinite number of distributions of current sources which can generate the same scalp electric potential and the same magnetic field outside the head and thus the inverse problem is highly ill-posed. This in turn necessitates the use of various assumptions which are often inappropriate; the technique also suffers from ‘cross-talk’ in that the activity at an arbitrary location is contaminated by activity from other brain regions.
One of the inventors finally provided a solution to the non-uniqueness problem, which is fundamental to reconstructing a current distribution based on MEG/EEG: although it is not possible to uniquely reconstruct the current in a conductor from knowledge of the electric potential on its boundary, or of magnetic field outside the conductor it is possible to determine which part of the current can be reconstructed via MEG and/or EEG. This solution is described in “Electro-magneto-encephalography for a Three-shell Model: Distributed Current in Arbitrary, Spherical and Ellipsoidal Geometries”, A S Fokas, J. R. Soc. Interface, 6:479-488, 2009. Further background prior art can be found in A S Fokas et al., “Electro-magneto-encephalography for the Three-shell Model: numerical implementation via splines for distributed current in spherical geometries”, Inverse Problems, 28, 2012. General background can be found in Liang et al., “Source reconstruction of brain electromagnetic fields”, Neuroimage, 27:1301-1311, 2009; Huang et al., “A novel integrated MEG and EEG analysis method for dipolar sources”, Neuroimage, 37:731-748, 2007; EP1468647; and JP2005/287675A.
The primary neuronal current distribution JP (T) is specified by three components of the current. This can be modelled by a combination of a scalar function, ψ, and a vector function A (altogether) with the constraint that the divergence of A is zero (giving 3 independent components). The use of EEG provides information about ψ, more particular the electric potential on the scalp can be represented by a combination of ψ and a second function vs (where s denotes the surface of the scalp). The use of MEG provides information about A; more particularly the MEG signal can be modelled by a combination of ψ and A, in combination with an auxiliary (vector) function H. Both H and vs are independent of the current distribution and depend only on the geometry and conductivity of the head.
The function H, is derivable from vs. In addition Fokas identified a set of partial differential equations whose solution defines vs (equations 2.11-2.14 in the paper, ibid). These are adumbrated below (renumbered as 0.11 et seq):
                                                        ∇              2                        ⁢                                          u                c                            ⁡                              (                                  r                  ,                  τ                                )                                              =          0                ,                  r          ∈                      Ω            c                          ,                  τ          ∈                      Ω            c                          ,                                  ⁢                                                            ∂                                                                                              ∂                n                                      ⁡                          [                                                1                                                                                r                      -                      τ                                                                                          +                                                      v                    c                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                              ]                                =                                    σ              f                        ⁢                                          ∂                                                      v                    f                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                                              ∂                n                                                    ,                              r            ∈                          ∂                              Ω                c                                              ;                                    (        0.11        )                                                                    ∇              2                        ⁢                                          v                f                            ⁡                              (                                  r                  ,                  τ                                )                                              =          0                ,                  r          ∈                      Ω            f                          ,                  τ          ∈                      Ω            c                          ,                                  ⁢                                            v              f                        ⁡                          (                              r                ,                τ                            )                                =                                    1                              σ                c                                      ⁡                          [                                                1                                                                                r                      -                      τ                                                                                          +                                                      v                    c                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                              ]                                      ,                  r          ∈                      ∂                          Ω              c                                      ,                                  ⁢                                            σ              f                        ⁢                                          ∂                                                      v                    f                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                                              ∂                n                                              =                                    σ              b                        ⁢                                          ∂                                                      v                    b                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                                              ∂                n                                                    ,                  r          ∈                      ∂                          Ω              f                                      ,                              τ            ∈                          Ω              c                                ;                                    (        0.12        )                                                                    ∇              2                        ⁢                                          v                b                            ⁡                              (                                  r                  ,                  τ                                )                                              =          0                ,                  r          ∈                      Ω            b                          ,                  τ          ∈                      Ω            c                          ,                                  ⁢                                            v              b                        ⁡                          (                              r                ,                τ                            )                                =                                    v              f                        ⁡                          (                              r                ,                τ                            )                                      ,                  r          ∈                      ∂                          Ω              f                                      ,                                  ⁢                                            σ              b                        ⁢                                          ∂                                                      v                    b                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                                              ∂                n                                              =                                    σ              s                        ⁢                                          ∂                                                      v                    s                                    ⁡                                      (                                          r                      ,                      τ                                        )                                                                              ∂                n                                                    ,                  r          ∈                      ∂                          Ω              b                                      ,                              τ            ∈                          Ω              c                                ;                                    (        0.13        )                                                                    ∇              2                        ⁢                                          v                s                            ⁡                              (                                  r                  ,                  τ                                )                                              =          0                ,                  r          ∈                      Ω            s                          ,                  τ          ∈                      Ω            c                          ,                                  ⁢                                            v              s                        ⁡                          (                              r                ,                τ                            )                                =                                    v              b                        ⁡                          (                              r                ,                τ                            )                                      ,                  r          ∈                      ∂                          Ω              b                                      ,                                  ⁢                                            ∂                                                v                  s                                ⁡                                  (                                      r                    ,                    τ                                    )                                                                    ∂              n                                =          0                ,                  r          ∈                                    ∂                              Ω                s                                      .                                              (        0.14        )            
In the above equations Ω denotes a three-dimensional space occupied by a conducting medium and ∂Ω its boundary. The head is modelled as four regions, a central region Ωc occupied by the cerebrum/brain (subscript c), around which are located three shells, Ωƒ, Ωb, and Ωs, respectively, modelling the spaces occupied by the cerebrospinal fluid (subscript f), the skull (subscript b), and the skin/scalp (subscript s). The region external to the head is denoted Ωe. In the equations ∂Ωc, ∂Ωƒ, ∂Ωb, and ∂Ωs denote the smooth boundaries of the domains Ωc, Ωƒ, Ωb, and Ωs, respectively. Each region has a respective electrical conductivity Υc, σƒ, σb, σs; the magnetic permeability is μ. The r vector (r∈R3) denotes position; later the (primary) neuronal current, the magnetic field, and the electric potential are denoted JP (r), B(r) and u(r) respectively. Similarly τ denotes a position vector within Ωc (where the primary current, JP (r), is defined). The derivative
      ∂                      ∂    n  denotes differentiation along the direction of {circumflex over (n)} where the hat denotes a unit vector and {circumflex over (n)} denotes the unit vector outward normal to the surface ∂Ω. As stated above, equations (0.11)-(0.14) are independent of the current JP(τ) and depend only on the geometry and on the conductivities σc, σb, σƒ and σs.
The skilled person will appreciate that to solve a set of partial differential equations the information needed comprises the equations, the respective domains of the equations, and the boundary conditions, as given above. Solving the above equations would lead to a value for vs, which could then be used to solve the inverse problem for EEG and/or MEG. However there is a problem in that the above equations cannot readily be solved.